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Abstract: In this paper we study the kinetic theory of many-particle astrophysical 
systems and we present a consistent version of the collisionless Boltzmann equation in the 
1PN approximation. We argue that the equation presented by Rezania and Sobouti in A&A 
354 1110 (2000) is not the correct expression to describe the evolution of a collisionless self- 
gravitating gas. One of the reasons that account for the previous statement is that the energy 
of a free-falling test particle, obeying the 1PN equations of motion for static gravitational 
fields, is not a static solution of the mentioned equation. The same statement holds for the 
angular momentum, in the case of spherical systems. We provide the necessary corrections 
and obtain an equation that is consistent with the corresponding equations of motion and 
the 1PN conserved quantities. We suggest some potential relevance for the study of high 
density astrophysical systems and as an application we construct the corrected version of 
the post-Newtonian polytropes. 
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1 Introduction 



The evolution of self-gravitating systems is a subject of great interest in astrophysics. From 
the statistical standpoint, the most straightforward way to construct self-consistent stellar 
systems is by means of finding the distribution function (DF) for a stellar system with 
a known gravitational potential and matter distribution. Since the mass density is the 
integration of the distribution function over the velocity variable in the phase space of the 
system, the problem of finding a DF is that of solving an integral equation (see [TJ [2j [3j HI [5] 
and the references therein). This construction is also the so-called "from p to /" approach 
for finding a self-consistent distribution function / [Bj, although the opposite procedure is 
also used sometimes. 

Now, in the framework of the general theory of relativity it is assumed that the DF 
satisfies the general relativistic version of the Fokker-Planck equation El E] or the colli- 
sionless Boltzmann equation (CBE) [10l[TT]. The former is devoted to systems in which local 
gravitational encounters play an essential role in their evolution whereas the latter is useful 
to study systems sufficiently smooth, so that they may be considered to be collisionless [6j. 
One can actually consider systems in which a number of particle species can collide and 
produce different species. This is how the formation of the light elements in the big bang 
nucleosynthesis is calculated (see [12] for a review). 

However, if we want to study the dynamics of huge astrophysical ensembles such as 
galaxies and galaxy clusters, physical collisions between the stars are very rare, and the 
effect of gravitational collisions can be neglected for times far longer than the age of the 
universe. Those systems are characterized by a relaxation time, t re i a x; that is arbitrarily large 
in comparison with their crossing time, t C ross, and this means that they can be approximated 
as a continuum rather than concentrated into nearly point-like stars. The same holds (with 
some restrictions) in the case of collisional systems such as globular clusters, neutron stars 
and withe dwarfs, where the relativistic effects of gravitation become important. Although 
^reiax here is significantly smaller than the system's age, the CBE is still valid over periods of 
time shorter than t re iax or when it is recognized that the system evolves slowly towards the 
equilibrium (on a timescale of the order of £ re iax)- For example, Taruya & Sakagami showed 
in [HI [13] , that the evolution of spherically symmetric systems in the collisional regime can 
be modeled as a sequence of polytropic states (i.e. described by a DF proportional to E 1 , 
which is a static solution of the CBE), with increasing polytropic index. 

There are many systems in astrophysics where Newtonian gravity is dominant, but general 
relativistic gravity plays also an important role in their evolution. For such systems it would 
be nice to have an approximation scheme which gives a Newtonian description in the lowest 
order and general relativistic effects as higher order perturbations. The post-Newtonian 
approximation is perfectly suited for this purpose. The appropriate scheme that describes 
the effects of the first post-Newtonian corrections on the motion of test particles, was first 
formulated by Einstein and Infeld [15} IT6j [T7] , then revisited by Weinberg [18] and it is known 
as the first post-Newtonian (1PN) approximation. This approach gives the corrections up to 
order v 2 / c 2 , where v is a typical velocity in the system and c is the speed of light, and it holds 
for particles moving non relativistically (v <C c), as in the case of a star moving around a 
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typical galaxy. Currently, PN approximations to higher orders have been developed because 
of the increasing interest around kinematics and associated emission of gravitational waves by 
binary pulsars, neutron stars and black holes. The most promising candidates for detecting 
such profiles are LIGO, VIRGO and GEO600 (see pH [20] for references). 

The first attempt to derive a kinetic equation in the 1PN approximation was done a few 
years ago by Rezania and Sobouti [21J, finding some relevant solutions. They tried to obtain 
the post-Newtonian version of the Liouville's equation for a gas of identical collisionless par- 
ticles and, as an application, they constructed the 1PN version of the classical polytropes PQ. 
Strictly speaking, this equation should be called the collisionless Boltzmann equation since 
the one attributed to Liouville refers to an iV-particle equation. However, both equations 
can be shown to coincide in the case N ^> 1 and when the iV-body distribution function is 
separable (see for example (6]). 

The authors started from the general relativistic kinetic equation 

£uf = (u^ - r^ir^-r) f{z*, u i ) = o, (i) 



where (x M , U l ) is the set of configuration and four- velocity coordinates, F 1 ^ are the Christof- 
fel's symbols, £u is the Liouville's operator and /(x M , U l ) is the one particle DF. Then, using 
the fact that the four-velocity of the particle satisfies the relation 

g^U'lT = ~c\ (2) 

they performed an expansion of Cjj up to the order (v/c) 2 , where v is the typical Newtonian 
speed. The resulting post-Newtonian approximation of ([1]) was: 

dt dx i dx i dv i 
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df 
— = 0, 



where <ft is the Newtonian gravitational potential, i/j and are the post-Newtonian potentials 
and / is now interpreted as a DF depending on the spatial coordinates x, the Newtonian 
velocity v, and the time t. 

As it was pointed out by the authors, one can verify that the above equation admits the 
static solution (dtp/dt = dtp/dt = and & = 0) 

E* = V + 0+(20 2 + ^)/c 2 , (4) 

which could be interpreted as the 1PN generalization of the classical energy. Moreover, for 
the case of spherically symmetric systems, there appear three additional integrals of equation 
©: 

l*= e ijk x^v k {l-<P/c 2 ), (5) 

which could also be considered as post-Newtonian extensions of the classical angular mo- 
menta. Thus, DFs depending on the integrals (JIJ and (111) would represent the 1PN statistical 
description of spherically symmetric systems in equilibrium. 
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However, we can verify that the "energy" given by 
equations of motion for a static source of gravitation 



is not an integral of the 1PN 



~dt 



-V 



+ 



4v 



+ — (v V0)-— V0. 



(6) 



Moreover, if we assume spherical symmetry for the fields <ft an d ip, we can also verify that 
the "angular momentum" (jSJ) is not an integral of motion of the above equation. These 
statements imply that there is not any correspondence between the microscopic description 
of motion given by (0) and the macroscopic (statistical) description involved in Since the 
former has been demonstrated to be correct (see [T3J [HJJ [T7J [TSJ dH [2H [22] for references), 
we can establish that the equation (J3J) does not provide an appropriate statistical description 
for the evolution of a collisionless gas of self-gravitating particles. Then, it is necessary to 
derive a corrected version of this relation. Such version must admit as a static solution, the 
correct expression for the energy at 1PN order (see section [2]): 



E 
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3v 4 3v 2 
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and, in the spherically symmetric case, the 1PN generalization of the angular momentum: 



(7) 
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In section [3] we show a detailed derivation of the corrected version of the CBE in the 1PN 



approximation. The final result reads 
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We can verify that the energy, given by (jTj), is a static solution of (jUJ), as well as, the 1PN 
angular momentum (jHJ), for the spherically symmetric case. It is worth to point out that the 
corrections introduced above, should be taken into account in the analysis of [23] and [24J, 
since they used some of the results of [21j. 

The rest of the paper is organized as follows. In Section [2] we start by deriving the 
post-Newtonian integrals of motion corresponding to static and spherically symmetric grav- 
itational field. Then, in Section [31 we show a detailed obtention of (EJ), and we rewrite it 
in a number of ways fl38l) - fl40l) in order to shed some light on its physical meaning. Finally, 
in Section HJ we construct the corrected version of the post-Newtonian polytropes. To do 
so, we employ the "/ to p" approach. We start with a DF that is just a generalization of 
the classical one and then we integrate over the phase space in order to obtain the different 
components of the stress tensor. Then we use the field equations to recover the Newtonian 
and post-Newtonian potentials, and we use them to study the main properties of the models. 
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2 Integrals of Motion in the 1PN aproximation 

We start by showing a detailed derivation of the first integral of motion, i.e. the energy, for 
a free-falling test particle in the presence of static gravitational fields (d<fi/dt = dip/dt = 
and £j = 0), in the 1PN approximation. In order to do this we consider the Lagrangian 
corresponding to the free-falling motion, 

22 = g^lT, (10) 

where g^ v is the metric tensor. We shall adopt the following conventions: 

f)'T^ ft HP F)'T'^ 

x , = (ct , x% I p = * = ** =d . v ,, ( „) 

and 

t = §L = V* = (MVc), (12) 

where r denotes the particle's proper time, Greek indices range from to 3 and Latin indices 
range from 1 to 3. 

At 1PN order, the line element ds 2 = g ilv dx il dx v , can be written in terms of the Newtonian 
potential <ft and post-Newtonian potentials ip and & defined as in [18] 



ds 2 = - (l + ^ + 2 ^ 2 + ^ c 2 dt 2 + 2 ^|^J cdtdxl 5 .. dx i dx i, (13) 

Then, for the case of static configurations, the post-Newtonian Lagrangian fflUT) reduces to 
2 C =-(l + ^ + ^)^ + (l-^)^, (14) 

where the upper dot denotes derivation with respect to r. It is possible to verify that, upon 
the application of the variational principle, the corresponding Euler-Lagrange equations lead 
to the 1PN equations of motion showed in (JH]). Since this Lagrangian is time- independent, 
the energy 

E=- 8 4 (15) 
dt 

is the first integral of motion. Our goal is to compute this quantity up to order c~ 2 . By 
taking the derivative of the Lagrangian with respect to t we obtain then 



20 2(0 2 + VQ 



B =( 1 + # + ^T^J Ct (16) 

Now we have to write t in terms of the spatial coordinates. In order to do this, first remem- 
ber that i = U°/c. The quantity U° can be computed by means of equation ([2]) in 1PN 
approximation, 



1 + 2| + 2(^W)> + (l-^ OT W, (17) 



and after some calculations we can find, up to order c 



-4 



U° _ v 2 <p 3v 4 5v 2 2 V _ ■ 



Finally, by introducing the above relation in (1161) we obtain the energy at order c 



-2. 



„ v 2 , 3v 4 3v 2 2 V 

Note that we have suppressed the constant term c 2 , which is meaningless for our purposes. 
By taking the total derivative of ( fl9i) with respect to t, and using the 1PN equations of 
motion ([6]), we can indeed prove that E is a conserved quantity at order c~ 2 . It turns out 
that this quantity is also a static solution of the corrected CBE (J5J), which can be proved by 
simple inspection (see section HJ eq. (jUJ). 

Using the 1PN equations of motion we can also rewrite the Lagrangian (j!4p . up to a total 
derivative, as 

which in fact agrees with the Lagrangian presented in [18] . Now, if we assume spherical 
symmetry, the fields and ip will depend on the spatial coordinates only through r = 
(x 1 ) 2 + (x 2 ) 2 + (x 3 ) 2 and, in this case, the Lagrangian has two cyclic coordinates 6 and <p 
given by 

tan<£> = x 2 /x 1 , (21) 
cos# = x 3 [(x 1 ) 2 + (x 2 ) 2 + (x 3 ) 2 ]" 1/2 . (22) 

This fact enable us to find that the three quantities 



k = £ ijk x>v k 1-4 + ^1 > ( 23 ) 



30 

c 2 2c 2 y 

are also integrals of motion and to interpret them as the post-Newtonian generalization of 
the angular momenta. For axial symmetry, of course, it is straightforward to show that only 
the ^-component of the angular momentum is conserved. 

The quantities (fT9l) and (J23]) are not static solutions of the equation (j3J), obtained in [2Tj . 
but they will be integrals of the corrected 1PN CBE, as expected. The derivation of such 
equation will be the aim of the next section. 

3 Derivation of the CBE in the 1PN Approximation 

In order to obtain a post-Newtonian approximation of equation ([TJ), it is convenient to change 
from (x^, U l ) to (x 11 , v l ). In other words, we use a transformation of the form 

(x^U 1 ) ^ (x^v'ix^U 1 )), (24) 
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so that the distribution function becomes 



(25) 



The explicit dependence of v 1 = v l (x fl , U l ) can be inferred from (j2j) and the fact that U l 
U°v l /c (see ( fTTT) and (JT2l) ). This allows us to write 



(26) 



To simplify the calculation we will first use the variable V 1 = (1, v i /c) and then we will come 
back to v % . The partial derivatives of / transform as 



df_ 

dx^ 

dj_ 



df_ 

dx^ 

df_ 



df_ 
dV j 



dV j 
dx 11 



(27) 



where the subscript at the bottom of the various derivatives indicates the quantity that is 
constant when / is differentiated. According to the above relations, equation ([T]) is rewritten 

as 



where 



df_ 

dx^ 
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df_ 



dV j 
dxi 1 
dV j 
dU 1 



dV j 
dx^ 



df_ 

dVi 



dV j 
dU 1 



U j fdU° 



u° 2 \dx^ 

5 ij U j fdU° 



U° U° 2 KdU* 



= 0, 
(28) 

(29) 
(30) 



The terms (dU° /dx^) vi and (dU° / dU 1 ) xil can be obtained by differentiating equation 
The result is 



dU° 
dx^ 
dU° 
dU 1 



■v a v p , 



U° 2 dg a ^ 
2Q dx» 

-^(gi V° + g ik V k ) 



and 



where 



Q = U\g 00 + g ol V l ) 



(31) 
(32) 

(33) 
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Thus, the required derivatives in equation ( 128|) are 

= U -V^V«V^ (34) 



<9v j ' \ V j 

"TTlfl'oi + ftfc^ ); for z 7^ j, 



-hc 2 U°~ 2 + J2v k (9ok + g k iV 1 )) for* 



(35) 



(note that {dV^ / dx^) ui differs from the expression shown in [2T] by a sign). After some 
calculations (see Appendix |A} one can verify that, up to order c -2 , the Liouville's operator 
in equation f )28|) can be expressed as 

C v = C d + C pn 

dt dx l dx l dv l c 2 \ 2 / \dt dx 1 



c z 



Av v J —— — h 30 7— + 3v — 7— — V J 

ox 1 \ 2 / ox 1 at ox 1 at 







\dx 3 


dx l ) 



d_ 

dv l 
(36) 

where £ is the classical Liouville operator (the first three terms of the r.h.s.) and C pn is 
the corresponding post-Newtonian correction (all terms multiplied by 1/c 2 ). Then, equation 
([28]) now reads 

(£ d + /7")/(x,v,t) = 0. (37) 

Thus we conclude that the CBE in the 1PN approximation can be splitted into a Newtonian 
contribution and a post-Newtonian one, as was obtained by [21] . but now with the corrected 
version for the operator C pn . 

Similar to the classical case, the 1PN equation ( 136]) can be expressed in various ways 
(see Appendix |B]) , each of which is useful in different contexts. First, as a vanishing total 
derivative, 

| = 0, (38) 

meaning that the flow through phase space of the probability fluid, (as seen by an observer 
moving with the particle) is incompressible j6]. Second, in terms of Poisson brackets, 

f + {/,#} = 0, (39) 

where H is the 1PN Hamiltonian. Since all integrals of motion must commute with H, it 
implies that Jeans theorem [25] is also valid at 1PN order, i.e. that any static solution of 
the CBE depends only on the integrals of motion of the system, and that any function of 
the integrals yields a static solution of the CBE. And third, as a continuity equation, 

! + ^<^°- < 4o > 
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where the set w = (q, p) is an arbitrary system of canonical coordinates. This equation 
states that the probability is conserved in phase space and, upon the appropriate integration 
over the momentum space, leads to the conservation laws in configuration space, i.e. the 
conservation of the energy-momentum tensor at 1PN order. We will refrain from writing out 
these results here, since they are not particularly illuminating. 



4 Static Solutions of the Post-Newtonian CBE 

For systems in static equilibrium / does not depend explicitly on time and the post-Newtonian 
potential & vanishes. In consequence, equation (J5J) reduces to: 



v 2 <p\ t d ( 3v 2 30 \ <90 d 4vV d<p d I dip d 



1 + - ^ V 1 — - 1 + -— + -M TTT ~n 7 + 



/(x,v) = 0. 



2c 2 c 2 / dx l \ 2 c 2 c 2 / dx % dv l c 2 dx % dvi c 2 dx l dv l 

(41) 

Since the energy (eq. ©) is an integral of motion of the system, we say that any ergodic 
DF f(E) satisfies (l4"Tj) . due to Jeans theorem. 

On the other hand, the gravitational fields and ip are related to the matter-energy 
distribution through the Einstein equations which, in the 1PN approximation, can be written 
as [18] 

V 2 = AnG °T 00 , V V = 4ttGc 2 ( 2 T 00 + 2 T U ) . (42) 
Here, the convention used is 

so that the symbol N T fll/ refers the /i^-component, of order (v/c) , in the expansion of 
the energy-momentum tensor. In particular, °T 00 is the density of rest-mass, 2 T 00 is the 
nonrelativistic part of the energy density and 2 T U (summation over i) is the classical kinetic 
energy density. 

Now, in general, the energy-momentum tensor is related to the DF through the equation 

T^(x x ) = -J ^rf{x\ U l )^d 3 U, (44) 

so that relations (|4ip -(l42 p form a set of self-consistent equations. For practical purposes it 
is necessary to expand f|44|) at various orders in v/c so, in order to illustrate this idea, we are 
now going to focus in the same special case used in [21], i.e. the post-Newtonian spherical 
polytropes. Needless to say, in the future it would be also interesting to consider the case 
of anisotropic systems, or axially symmetric systems, specially to develop applications for 
galactic dynamics. 
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4.1 Construction of Post-Newtonian Polytropes 

Now we deal with the 1PN version of polytropic solutions, i.e. systems characterized with 
DFs of the form 

f(E)=^(-E)^ 2 for£<0, 

JK J 2n J ' (45) 

= for E > 0, 

where k n is a real constant and n is the index of the polytrope. The DF can be splitted in 
two parts: a c° (classical) contribution /(°) and a c~ 2 (post-Newtonian) contribution f( 2 \ 
This is possible because the energy can be written as 

E = E d + E pn (46) 

where 

v2 j. rn 3v 4 3v 2 (j) 2 tP 

Since we assume that E d 3> E pn , we can write 

/ « ^ (-^d)"" 372 [1 + (n - 3/2) ^ n /E d ] 

= ^ (-^)"~ 3/2 - ^ (n - 3/2) £ pn (-^)"" 5/2 ^) 
= /°> + / (2) 

In order to obtain T M1/ at various orders of v/c, we will need the expansion of U° (see eq. 
(fl8|) ). the expansion of the determinant of the metric tensor, y/—g = 1 — 20/c 2 + ... and also 
the expansion of d 3 U, rewritten in terms of d 3 v. Remembering that U % = U°v l /c and using 
(ITS]) , we have 



dU* d 





v 2 


(- 


2c 2 


v 2 


<\> 


2^ " 


c 2 


v 2 


<\> 


2c 2 " 


c 2 



c 2 



1^-51^4(^1 (49) 



<9t> J V 2c 2 



which means that the Jacobian of the transformation is \dU l /dv^\ = l+5v 2 /(2c 2 ) — 30/c 2 + ... 
and, in consequence, 

<ft7=(i + g-£)A. = *r(i + g-£)v* (50) 

(the r.h.s. of the last equation can be implemented only in the case of a DF depending on 
the velocity components through v = (fjf*) 1//2 ). Putting all of this together, we can now 
write T 00 and T n as: 

T 00 = 4vr £ [/(•») + (l + ^ " ^) v 2 rf. (51) 
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and 



I*=£jf [/<°»+/ (2 '](l + ^-f)v^. (52) 

The components T°* = T'° vanish due to the distribution of matter is static (see [H]). Here 
v e denotes the escape velocity, i.e. the speed at which a particle reaches its maximum value 
of energy, E = 0, so that it is confined to the distribution of matter. Such quantity can be 
computed from ([7]) and the result is 



, 2c 2 2 
v e =\l + 20 + - x /c 4 -12c 2 + 60 2 - 6?/> 
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(53) 



In order to be consistent with the 1PN approximation, one can set v e = y/— 20 due to the 
fact that the c-dependent terms would lead to the apparition of a factor of order 0(c~ 6 ), 
even in the case when the integrand is f(°\ Introducing this value in ()5TT)-( l52l) . we can obtain 
explicitly °T 00 , 2 T 00 and 2 T U (see Appendix [C]) and, after some calculations, we obtain 

V 2 = a„(-0) n , (54) 
V 2 ^ = -na n (-<P) n - 1 il> + P n (-<f>) n+1 , (55) 

where V 2 = (l/r 2 )(d/dr)(r 2 d/dr) and we have introduced the constants 

a n = W2 ^Gk n r l n .-y^\ (56) 

r(n + 1) 

r w>„, (n 2 - 2n - 63/4)r(n - 1/2) 

(3 n = -2V2 ^' 2 Gk n y - 1 / y '-J-. 

lin + 2) 

The first equation above fl54l) (which looks like very different from the one derived in (21]) is 
the classical field equation for the Newtonian polytropes [lj and has simple exact solutions 
for the cases n — 0, 1, 5, the latter corresponding to the Plummer's model [2S1EZ1- For other 
values of n, the solution can not be expressed in terms of elementary functions [6] . 

4.2 Numerical Solutions of the Field Equations 

In order to perform the numerical solution of the system for any n, we implement the 
following definitions: 

~ / — ? j. y 1 n (n 2 - 2n - 63/4)02 

r =\/a n (-6 ) n - 1 r, B n = ; . / >Yo , 

V 1 V ' ' 2(n + l)^ (57) 

=0 O X, ip = ip Y, 
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where <p Q and ip Q are the Newtonian and post-newtonian gravitational potentials at the center 
of the configuration, respectively. Then, relations (I54l) - (j55j) become 

X n , (58) 
nX n ~ x Y - B n X n+1 . (59) 

Since we assume that the gravitational potentials reach critical values at the center of the 
configuration, we have to impose the initial conditions 

X(0) = Y(0) = 1, X'(0) = Y'(0) = 0, (60) 

where the prime denotes differentiation with respect to the scaling radius f. We use a fourth- 
order Runge-Kutta method to find the numerical solutions for different values of the ratios 
<po/c 2 and assuming ip D ~ — <p 2 . In a neutron star, for example, <p /(? can vary between — 1 
and —0.1, approximately. 

In Figure [1] we plot the gravitational potential energy of a test particle with velocity 
v — (scaled with respect to <p a ), given by 

for the case of Newtonian polytropes and post-Newtonian polytropes. We note that for 
<p D /c 2 ~ —0.02 or smaller, the post-Newtonian corrections are not significant, while for 
values <p /c 2 ~ —0.2 or larger they become important. 



1 


d , 


> „ 2 dX 
\ df 




df 


1 








df 



4.3 Post-Newtonian corrections to the rotation curves and mass 
densities 

In this subsection we study the fundamental equations describing the circular motion of test 
particles, in order to investigate the corrections introduced by the relativistic effects on the 
rotation curves. In spherical coordinates (r,8,ip), the circular orbits in the equatorial plane 
can be obtained from (EJ), by setting 9 = tt/2, d(p/dt = and dip /dip = dip /dip = 0. They 
must satisfy the conditions f — z — 0, r—'z — O and z = 0, so that the equation of motion 
reduces to 



rip 



~^ c 2 dr 



d_ 

dr 



+ 



24> 2 + iJj 



This can be used to derive an expression for the circular velocity v v 
the radius r. The result reads: 



40 r dd\ r dip 
1 + -it o ~x + -^-^r- 



dr 



c 2 dr 



c 2 dr 



(62) 

rip, as a function of 
(63) 



z=0 



Note that in the limit c — > oo, the above expression reduces to the usual relation derived in 
Newtonian theory: v v = \J Rd<p/dR. Perhaps, the most important difference between such 
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relation and ( )63|) is that, in the Newtonian case, the radical is linear in and its derivatives, 
whereas in the 1PN case, it depends on non linear terms involving 0, if> and derivatives. This 
non linear dependence may result significant in some cases and its effects can be observed 
in the rotation curves. In Figure [2] we show the rotation curves corresponding to the 1PN 
corrected models. As we can see, in some cases the 1PN corrections are significant while, in 
other cases, they are practically negligible. 

Now, in order to examine if the 1PN corrections provide an adequate physical description 
for the energy-mass distribution, we plot the density p. The contribution of °T 00 plus 2 T 00 
(see Appendix [C] for details) as a function of r is shown in Figure [3] (here p is the mass 
density divided by k n (f>l, in each case). We find that p is a positive- valued function with a 
maximum in the center and a minimum at r — > oo, for the following situations: (i) polytropes 
n = 1, . . . , 6 with (j) /c 2 = —0.02; (ii) polytrope n = 6 with (p /c 2 = —0.2. When such ratio 
reaches the value —0.4, the mass density becomes negative for certain values of r, far from 
the center (see Figure H]), thus representing a non-physical situation with tachyonic matter. 

5 Concluding Remarks 

We have obtained the 1PN version of the collisionless Boltzmann equation for a self-gravitating 
gas of identical particles, which is consistent with the microscopic equations of motion. This 
can be shown by checking that the integrals of motion derived from the microscopic dynam- 
ics, are static solutions of (Q, in agreement with the statistical description of the system. 
Such relation leads to the equations of the post-Newtonian hydrodynamics, derived previ- 
ously by Chandrasekar [22], and implies automatically that the macroscopic stress-energy 
tensor is conserved if one takes the corresponding integrals over the phase space. 

The interpretation of the CBE derived here follows the same logical arguments as in the 
classical case: (i) we can think of it as a continuity equation of the probability fluid, which 
is just a statement that comes from the probability conservation in phase space, (ii) The 
fact that we are allowed to rewrite it as a vanishing total derivative means that the flow of 
this probability fluid, as seen by a comoving observer, is incompressible, (iii) The structure 
of the equation in terms of poisson brackets reveals that Jeans theorem is also valid at 1PN 
order, i.e. that any static solution of the CBE depends only on the integrals of motion of 
the system, and that any function of the integrals yields a static solution of the CBE. We 
suggest some potential relevance of our findings for the study of nuclear cores in galactic 
dynamics and other astrophysical systems with high enough energy densities such that the 
relativistic effects become important. 

As a first step towards the developing of such astrophysical applications, we considered 
the case of post-Newtonian polytropes, thus providing the corrected version of the solutions 
found in [21] . We proceeded to do it numerically because, even in the Newtonian case, it is not 
possible to find analytical solutions for arbitrary polytropic index. Now, in order to compare 
the behavior of Newtonian and 1PN solutions, we chose some values for the dimensionless 
parameters 4> /c 2 and ip/4>l, so that they can be associated to relativistic stellar systems. 
For example, for (p /c 2 = —0.02 and ip/4> 2 = —1, we found that the Newtonian and post- 
Newtonian behavior are very similar, although some perceptible differences can be seen in 
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the gravitational energy, rotation curves and mass densities. In particular, we note that 
that the 1PN corrections to the mass densities are more significant near the center of the 
configurations and, in contrast, the 1PN corrections to the circular velocities are greater far 
from the center. In general, the 1PN values of v c are smaller that the Newtonian values. 
When we chose 4> /c 2 = —0.2 or 4> /c 2 = —0.4, we found that the relativistic contributions are 
more relevant. However, in some cases the 1PN models are unphysical for these parameters 
(meaning that the mass density becomes negative for some values of r), depending on the 
polytropic index n. We also note the fact that, the larger are the corrections to the central 
mass density, the larger are the corrections to the circular velocity far from the center. 

Last but not least, it is worth to point out that it is also possible to derive an exact 
solution for the case n = 5, i.e. the 1PN version of Plummer's model, which is the simplest 
polytrope with physical relevance. It will be presented in a subsequent paper with further 
applications. 

Acknowledgements 

J.R.-C. is grateful to FAPESP for financial support and to Prof. P. S. Letelier for stimulating 
discussions. 



A Derivation of relation (1361) 



By introducing (|34|) and (I35I) in the equation (128|) . the Liouville's operator takes the form 

Cu ~ UV + ~2Q-^ V VV WJ- -Q~ r ^ V V ^ -Q i9oi + 9ikV ] dVi 





(U ) 2 ■ 

+ ^ — ! —r i vv 1 

Q » v 



(U°r 2 + J2 vk (9 k + 9kiV l ) 



dV i 



(64) 



The four terms in the r.h.s., up to order c 2 , can be written as 

wJ- = (l + *-+)° + (l + *-+) v *±, 
8x» V 2c 2 c 2 dt \ 2c 2 c 2 dx v 



(65) 



v l d4> d v l vi d(t> d 



2Q dx^ dVi c 2 dt dv l c 2 dx l dvi 



(66) 



P°) 2 ri/T V F fo + a V k ) 9 9 {VI)29(P 9 (67) 
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Q 



c 2 (u°y 



3v 2 



J2 Vk (9ok+9klV l 



1 + 7^-7 + ^ 



2 c 2 
2v l d(f) d 
c 2 dt dv 



i\2 



3<A d(f) d (v l ) 



2 / dx l dv 1 
1 dip d 
c 2 dx l dv 1 



+ 



d_ 

J dV 1 





c 2 dx % dv % 



l_d^d_ 

c 2 dt dv 1 



dii da\ d 



c 2 V dxi dx l J dv 



In the derivation of (l65"]) - (l6"8j) we have used the approximations 



1 - 



+ 



(U 



0\2 



cQ 
c 2 (U°) 



•1 + 



3<f) 



2c 2 ' 

2^' 



1 + ^-^7- 



2v l v^ d(f) d 
c 2 dx l dvi 



(68) 



(69) 



Finally, replacing ( 165]) - (168]) in (164"]) we obtain the final form of the Liouville's operator, eq. 



B Alternative expressions for the CBE 

In order to express the CBE as a total derivative, let's start by writing explicitly the various 
orders of the DF, i.e. 

C v f = C cl f w + (£ c 7 {2) + C pn f {0) ) = 0, (70) 

where 

£df(o) = o (71) 

and 

£dj(2) + £jmy(0) = q ( 72 j 

The first of these equations implies that 

d f(o) 9f (0) d(f)dfi o) 

+ v "17777 = 7^7^77- • ( 73 ) 



Moreover, 



whereas 



dt dx % dx l dv 1 



£ ^,^ + ^_M (74) 

dt dx % dx % dv 1 



i 

+ 



c 2 



4wV-^- - I — + 30 ) — ! °>r'— - — - r> 

dx ] V 2 / dx % dt dx % dt 
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\dxi 


dx i J 



dv 1 



Replacing (173]) in the previous equation we find 



1 



4-uV- 



-— - (v z + 40)-— + 3f 

ftr^ v 7 dx l 



dtp dtp d^i 
~dt ~ dot ~ ~dt 



dt ± _di 1 

dxi dx l 



df 



(o) 



Here we need the 1PN equations of motion which, in general, are given by 



(75) 



1 



4t> V - 



dx l c 2 

Introducing (1751) in (1751) we obtain 

£pn j?(0) 



v 2 + 40)^ + 3„i^ ^ g& t j 
cte* 1 dt dx l dt 



dxi dx l 



dt dx i 



df<® 
dv l 



(76) 



(77) 



so 



dt ' " dx l ' dt du* ' ^ 
Also, note that using the 1PN equations of motion in (174j) . we only have to keep the leading 



term: 



K2)_d/ (2) | vl dfM | dv'df 
dt dx % dt dv l 



(2) 



(79) 



Putting all of this together we can finally write 

£ v f = C cl f {0) + (C cl f {2) + CTf®) 



(80) 



or 



I- - < 81 > 

Now, if we use a set of canonical coordinates w = (q, p) (i.e. a set of coordinates that 
satisfies the Hamilton equations), we can rewrite ( 1HT1) as 



or 



which means that 



# = 3/ d_i_df_ di_df_ = dj_ 

dt dt dt dq l dt dt iJ ' 1 



di = 21 d ( m\ of df 

dt dt dq { V dt J dp* V dt 



df d 

^7 + 



dt dw r 



f 



dw n 
~dt 



(82) 
(83) 
(84) 



where m = 1, 6 and w m is the m-th component of w. 
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C Derivation of Equations (1541)- {55} 



In order to calculate the r.h.s of equations (I42p . we have to write the required components 
of the energy-momentum tensor, °T 00 , 2 T 00 and 2 T n (sum over % — 1,2, 3) explicitly 

r V=2? / 2 \ n-3/2 

°T 00 = 2*^ (-y-^J (85) 

fi r, fyf=2$ / 2 \ "-3/2 

2 t- = ^jf (-y-^) (v 2 -20)v 2 ^ 



20 / Ar 2 \ "-3/2 



2t " - %L ' {-T-+) ^ (86) 

The above formulae can be computed easily by making the substitution v 2 = —20 cos 2 9 and 
running the integrals over 9, from to 7r/2. The resulting trigonometric integrals can be 
expressed in terms of gamma functions by using the relation [28] 

*V 0)-(sin fff^O = r{k + m r + )'?, (87) 
,, J K J 2r(m + k + 3/2) ' v ; 

for m > —1/2 and k > —1. Then we obtain, for n > 1/2, 

°T oo = 2£; n /(0,n- 3/2,1), 

2 T oo + 2 T u = _n ^ n _ 3 ^ ^ _ 12 0j(0 ) n _ 3/2, 1) + (n - 3/2) [30/(0, n - 5/2, 2) 



-(0 2 + 2^)7(0, n - 5/2, 1) - (3/4)7(0, n - 5/2, 3)] 



where we have introduced the notation 



w ' ' r(m + k + 3/2) v r; ' v ; 



and after some calculations we obtain the expressions 
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Figure 1: We show the gravitational potential energy U (given by eq.( )6Tl) ) for Newtonian 
(continuous line) and Post-Newtonian polytropes with index n = 1, .., 6. In each illustration 
we assume ip a ~ 4> 2 and choose (p /c 2 equal to —0.02 (dashed line), —0.2 (dotted line) and 
-0.4 (dash-dotted line). 
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Figure 2: We show the circular velocity V c for Newtonian (continuous line) and Post- 
Newtonian polytropes with index n = 1,..,6. In each illustration we assume ip ~ <p 2 
and choose 4> /c 2 equal to —0.02 (dashed line), —0.2 (dotted line) and —0.4 (dash-dotted 
line). 
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Figure 3: We plot the mass density p, scaled with respect k n 4>®, for Newtonian (continuous 
line) and Post-Newtonian polytropes with index n — 1, ..,6. In each illustration we assume 
ip ~ (f) 2 and choose o /c 2 equal to —0.02 (dashed line), —0.2 (dotted line) and —0.4 (dash- 
dotted line). 
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Figure 4: Details of figure |3]for Newtonian (continuous line) and Post-Newtonian polytropes 
with index n = 1, ..,6. In each illustration we assume ip Q ~ <p 2 and choose <p /c 2 equal to 
-0.02 (dashed line), -0.2 (dotted line) and -0.4 (dash-dotted line). 
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